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The hydrogen plasma is studied at temperatures T ~ 
10* - 10^ K using the free energy minimization method. A 
simple analytic free energy model is proposed which is accu- 
rate at densities p ^ 1 g cm~^ and yields convergent inter- 
nal partition function of atoms. The occupation probability 
formalism is modified for solving the ionization equilibrium 
problem. The ionization degree and equation of state are cal- 
culated and compared with the results of other models. 

PACS numbers: 95.30. Q, 95.30.T, 52.25.K 



I. INTRODUCTION 



at high density. This compelled the authors to introduce 



compi 

ad hoc modifications |^,|lO| which affected the ionization 
equilibrium. The ionization curves obtained in Refs. |^ 
and at T > 10^ K are strikingly different. 

In this paper we study hydrogen plasma at T ~ 
10'* — 10^ K when a considerable fraction of atoms can 
exist in excited states, and the IPF convergence becomes 
crucial. In Sect. || we outline the free energy minimiza- 
tion method. In Sect. [II we develop an analytic free 
energy model for weakly coupled and weakly degenerate 
plasma. In Sect. ^ we revise the occupation probabil- 
ity formalism, and in Sect. ^ we calculate the EOS and 
ionization equilibrium, making a comparison with the re- 
sults of other models. 
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Thermodynamic properties of astrophysical plasmas 
have been studied extensively in recent years [Q. The 
theoretical models are based either on the physical pic- 
ture or on the chemical picture of the plasma . In this 
paper we consider the latter. Our study is motivated by 
considerable disagreement among the existing models in 
the domain of partial ionization. 

We consider the simplest case of pure hydrogen plasma 
which shows in relief all specific features of the prob- 
lem. This particular case is also important for calculat- 
ing atmospheric opacities of degenerate stars. It is gener- 
ally assumed that, due to the gravitational stratification 

, outer spectra- forming layers of atmospheres of these 
stars consist of light elements. Observations Q confirm 
that DA white dwarfs have virtually pure hydrogen at- 
mospheres at temperatures up to 40 000 K. 

It is worthwhile to mention two widely used theories 
for the hydrogen plasma. The first one has been elabo- 
rated by Hummer and Mihalas Q (HM), and Mihalas, 
Dappen, and Hummer |^J^ (MDH). The authors pre- 
sented a detailed discussion of previous work and formu- 
lated an equation of state (EOS), using an occupation 
probability formalism to obtain a finite internal parti- 
tion function (IPF) of atoms. The occupation probabil- 
ities were derived from the plasma microfield distribu- 
tion. The second theory has been proposed by Saumon 
and Chabrier (SC) (see also Ref. 0). These au- 

thors have developed a free energy model based on ef- 
fective pair potentials in the system of hydrogen atoms, 
molecules, protons, and electrons. The theory describes 
successfully either a completely ionized plasma or neu- 
tral gas 1^ but suffers from difficulties in treating the par- 
tial ionization. The free energy models originally adopted 
by HM 1^ and SC |^ did not ensure reasonable conver- 
gence of IPF at high temperature and pressure ionization 



II. CHEMICAL PICTURE AND FREE ENERGY 
MINIMIZATION METHOD 

In the chemical picture of plasmas, bound objects 
(atoms, molecules, ions) are treated as elementary mem- 
bers of the thermodynamic ensemble, along with free 
electrons and nuclei. In the physical picture, nuclei and 
electrons (free and bound) are the only constituents of 
the ensemble. 

Both pictures can be thermodynamically self-con- 
sistent, but the chemical picture has limited microscopic 
consistency. For instance, it does not provide a proper 
treatment of such cluster configurations as "an atom -I- a 
close alien ion" (a pseudomolecular ion) . With increasing 
density, the ionization of an electron bound to a particu- 
lar nucleus proceeds through a progressive delocalization 
involving cluster ("hopping") states [Q. These states 
are negligible in nearly ideal gases, but important at 
higher densities. However, their inclusion as new mem- 
bers of the thermodynamic ensemble would complicate 
the free energy model. Therefore one usually considers 
basic chemical species which dominate at low density. 

On the other hand, the physical picture is commonly 
based on diagrammatic expansions which converge only 
at low densities. Thus one has to resort to additional 
assumptions in the frames of either picture, in order to 
progress on higher density. The chemical picture, com- 
bined with the free energy minimization method, rep- 
resents a reasonable compromise between the rigorous 
treatment and the practical application. 

The central assumption of the free energy minimiza- 
tion method is the factorization of the many-body 
partition function into translational, configurational, and 
internal factors, and corresponding separation of the 
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Helmholtz free energy F: 



C. Internal free energy 



F{V,T,{Na,})^Fu 



^conf ~t~ ^ir 



(1) 



Here V is the volume, T the temperature, and {Na} the 
set of particle numbers. The internal structure of a com- 
posite particle is generally affected by the surrounding, 
hence the separation (|l|) is approximate. 

At given V and T, the equilibrium state is determined 
by minimizing F with respect to the numbers {Na}, sub- 
ject to the stoichiometric constraints. Then the pressure 
P, the entropy S and related quantities are obtained from 
F using the well known thermodynamic relations JT^ . 



III. FREE ENERGY MODEL 
A. Plasma parameters 

Consider a plasma of electrons, protons, and hydrogen 
atoms. We do not include molecules and molecular ions, 
assuming that the temperature is high enough for their 
dissociation (roughly, T > 10^ K; possible departure from 
this assumption will be discussed in Sect.^. 

The charged component of the plasma is described by 
the coupling parameter T and degeneracy parameter 0, 



iSeya,, 9 = T/Tf, 



(2) 



where /3 — (kBT)^^, fee is the Boltzmann constant, 
flo = {^TTric)^^^^ is the mean inter-electron distance, 
ric ~ NcjV v& the electron number density, and Tp 
is the Fermi temperature (Tp 
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[12 500K/T]n2^'', where nao = ne/lO^o cm'^). We con- 
sider weakly coupled and weakly degenerate plasma. The 
generalization to higher degeneracy (0 < 1) is straightfor- 
ward (l4j. The strong Coulomb coupling (F > 1) can be 
taken into account using the models developed in Refs. 



B. Translational free energy 

The free energy of the ideal plasma F\a — Ptrans+Pint is 
the sum of three terms produced by electrons (e), protons 
(p), and H atoms. For the atoms and protons. 



/S^'^Va^p = ^iF'^L|N^. = ln(r^pA3) - 1, 



/JFtS.s/A^H = In(nHA^) - 1, 



(3) 
(4) 



where Ua — Na/V is the number density of species a, 
and Aq, = {2it[3Ti /toq)^^^ is the thermal de Broglie wave- 
length. We neglect the proton spin weight since it would 
yield an insignificant constant in the free energy. For the 
electrons, we include the spin weight and the low-density 
correction for the degeneracy u%: 



^F^^^/Ne = Hn,Xl/2) - 1 + neA3/2^/2. 



(5) 



There are different ways to define the internal free en- 
ergy of atoms f /iJ^'' . We calculate it in the ideal gas ap- 
proximation, neglecting interactions of atoms with sur- 
rounding particles. From the first principles [Ti| ], the 
ideal-gas part of the free energy is 



(6) 



where n enumerates quantum states with statistical 
weights and non-perturbed binding energies Xk- It has 
been shown JlGf that the binding energies of an atom in 
a plasma practically do not shift with increasing density 
until they merge into the continuum. Comparing Eqs. 
M) and (H), we obtain 



(7) 



Now all nonideality effects should be included in the con- 
figurational term Pconf- However, Fconf does affect the 
equilibrium value of F^^^ through the distribution of the 
occupation numbers {A^k}, which is not assumed to obey 
the ideal-gas Boltzmann law. 



D. Configurational free energy 

It is the common practice to separate Coulomb inter- 
actions of charged particles from interactions involving 
neutral atoms and to describe the first ones by the free 
energy of a pure Coulomb plasma. At low density, the 
excess free energy of the charged component is given by 
the Debye-Hiickel theory with a two-component plasma 
quantum correction [T^ : 



PFc 



2e^ 



V 



2V2 



(8) 



where 7 = jijcc + 27cp), 7cc and 7cp being the electron- 
electron and electron-proton quantum diffraction param- 
eters, 7ep « 7cg/V2 = 2hf3e ^7r(nc -I- np)/me. Our 7 
differs from the one-component plasma parameter 7e em- 
ployed in the SC model [|| by multiplier 0.854. 

The quantum corrections in Eqs. (|^) and (j^) allow us 
to extend the analytic free energy model from the low 
density region of F ^ 1 and 6* ^ 1 to the moderate 
density, where F < 1, 6* > 1. We have checked (see also 
Fig. 12 of Ref. ||^) that the account of these corrections 
extends the validity range of the model by more than 
an order of magnitude toward higher density or lower 
temperature. 

The neutral component produces additional configu- 
rational terms, which describe interactions of neutral 
species with neutral and charged particles. For atoms 
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and ions, these interactions are described by pair poten- 
tials, while the electrons are assumed to adjust to any 
configuration of the heavy particles. A free energy model 
based on this approach has been elaborated by SC [p|Jlo| 
and extended recently to helium plasma p8[ . Our modi- 
fied version of the pair potential free energy model -Fpair 
is given in Appendix. However, the contribution to the 
free energy, which comes from unbounding of atoms in 
the course of their interactions with other particles, is 
not fully griped at this approach. Namely, the effective 
hard sphere diameters of atoms, derived from effective 
binary potentials, appear to be too small to enforce ion- 
ization in high-T and high-p domain. Below we consider 
the problem in more detail. 



To allow for the plasma screening of the intra-atomic 
Coulomb interactions, the static screened Coulomb po- 
tential (SSCP) has been widely used (e.g., Ref. 
However this approach has been strongly criticized[y . 
More consistent is the dynamical screening recently in- 
vestigated in the framework of the thermodynamic Green 
function method ||T^. However, the practical need (e.g., 
in astrophysics) for large arrays of thermodynamic data 
necessitates looking for a simplified approach. 

A convenient model has been described by HM. The 
configurational term given by Eq. (4.72) of Ref. ^ for 
the hydrogen plasma may be represented as 



(9) 



E. Pressure destruction of atoms 

The central problem of calculating the occupation 
numbers {Nq.} is to achieve the self-consistent conver- 
gence of IPF, that is most difficult in the case of partial 
ionization. There have been many attempts to solve the 
problem; some of them are discussed in Ref. [D . 

Clearly, a bound state cannot be populated, if the cor- 
responding "size" of the electronic "orbit" (the electron 
cloud) is too large. In other words, the electron cannot 
be bound if there is not enough free space, or if it suf- 
fers from strong perturbations. The problem is how to 
include the unbounding into the chemical picture in a 
self-consistent manner. 

The simplest way to obtain the convergence is to 
truncate the IPF at an appropriate "critical" quantum 
number k*, for example, corresponding to a quantum- 
mechanical size of an atom 1^' comparable with the mean 
interparticle distance However, since k* depends 

on the physical conditions (particularly, on density), the 
abrupt cutoff produces unrealistic discontinuities in the 
free energy. Continuous truncation procedures imply in- 
troduction of occupation probabilities w,^ into the IPF, 
which suppress higher states and ensure the convergence. 
However, as has been shown by Fermi |20[| and empha- 
sized by HM 1^ , the introduction of w„ requires a modi- 
fication of the free energy. Therefore tu^ should be con- 
sistent with the adopted form of i\:onf • 

In the free energy model Fpair (Appendix), the re- 
pulsive interparticle interactions are simulated by repul- 
sion of hard spheres (HS) with appropriate diameters d. 
However, for a high T or for a high atomic level, d be- 
comes much smaller than the atomic size so that other 
plasma particles can penetrate the atom. If this happens, 
then the electronic orbit gets embedded into the plasma, 
which screens the attraction to the nucleus and makes the 
electron unbound. Then the atomic constituents should 
be extracted from the neutral component and treated 
as independent participants in the charged component. 
Similar situation (known as pressure ionization) occurs 
even for ground state atoms, if the density is high enough. 



where vik and Wp^ are characteristic volumes associated, 
respectively, with the atom-atom and atom-ion interac- 
tions. The first term includes the interaction of an atom 
in a state k with the ground state atom only. This is the 
"low excitation approximation" proposed by HM to make 
^conf linear in , which was essentially employed in de- 
riving the occupation probabilities (despite the fact that 
even at this approximation the linearity obviously breaks 
down for the ground state atoms). For purely neutral gas 
(rip = 0), Eq. (^) reproduces the free energy derived by 
Fermi |20[| . The atom-ion interaction volume VpK, in this 
approach is due to microfield perturbations; it depends 
on the principal quantum number in a complicated way, 
but it is always 1-2 orders of magnitude larger than vik- 
For example, for the ground-state hydrogen, Eq. (4.70) of 
Ref. II yields Vpi = 128 (47r/3)a|, where cb is the Bohr 
radius. 

It was found Q| that this model does not provide 
pressure ionization at high density, which is not surpris- 
ing. Indeed, since Vp^ 3> wik, Pconf"^ may increase more 
rapidly with growing Hp than with njj, shifting the equi- 
librium towards lower ionization degree at higher densi- 
ties. In order to ensure the desired ionization, MDH in- 
troduced an artificial "pressure ionization term" F5 into 
the free energy, which rapidly increased whenever the 
density of neutral fraction exceeded 10~^ g cm~^. It can 
be shown that it is this term (and not the inaccurate 
HS treatment, as supposed in Ref. ||ll|) that produces an 
unrealistically stiff EOS at p > 10^^ cm~'^. 

As argued in Ref s. |n^ , pO| ,pt, the quantum- mechanical 
atomic size [Eq. ([A9|y | should not exceed the mean in- 
terparticle distance. This is not a problem at low tem- 
perature, when the effective HS diameters d^K' are larger 
than /k- At high T, however, d^K' become small, allowing 
configurations with strongly overlapping wavefunctions 
of neutral atoms. SC escaped this difficulty by introduc- 
ing an additional hard core [Eq. (14) of Ref. |10|] in the 
effective potentials. However, a large hard core inside the 
atom seems to be unrealistic. We propose a modification 
of the free energy, which has another interpretation. 

An electron can be treated as bound to a particu- 
lar nucleus, if only its wave function does not overlap 
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strongly with wave functions of other electrons, either 
free or bound to neighboring atoms. From the classi- 
cal point of view, the atomic electron becomes unbound 
when another electron falls inside its orbit and shields 
the attraction to the nucleus. This resembles the plasma 
screening of the nucleus in the SSCP model but does not 
imply the collective nature of the screening. 

For a given state k, the probability that such un- 
bounding does not occur can be estimated at low den- 
sity from the Poisson distribution, p„ = exp[— 
where = (tt-o + hh) is the total number of randomly 
distributed electronic clouds (including those which are 
bound to nuclei), and = |7rZ^. The unbounding re- 
quires to exclude the overlapping configurations from the 
total partition function, thus reducing the volume of the 
phase space available to the system. Equivalently, the 
existence of an atom in the state n corresponds to an 
event with probability and thus diminishes the en- 
tropy. The total negentropy corresponding to a set of 
occupation numbers {N^,} is — J2k Ii^Pk, which gives 
the free energy contribution 



Fuh = kBTNcnuvu, 



(10) 



where ?;h = '^h^kVk/N-h is the average atomic volume. 

A similar term has been introduced by HM who, how- 
ever, considered the destruction of atoms by microfields 
fluctuating due to the motions of surrounding heavy par- 
ticles (as discussed in Sect. IV6(ii) of Ref . [^). Although 
we readily agree that nearby passages of positive ions can 
ionize a particular atom, this process does not affect the 
occupation numbers at the thermodynamic equilibrium, 
since it is compensated by the inverse (neglected) pro- 
cess, owing to the principle of detailed balance. The net 
effect of both processes is not given in advance but itself 
should be determined from the thermodynamic equilib- 
rium conditions. In contrast, the unbounding by an ex- 
cessive negative charge occurring inside an electron orbit 
has no balancing counterpart. Moreover, it seems incon- 
sistent to include any dynamical process dependent on 
particle momenta into -Fconf after separating the trans- 
lational term -Ftrans, because the separation (^ implies 
that the other terms may depend only on particle config- 
uration coordinates. For this reason we do not consider 
atomic collisions with ions Q and free electrons [Q. 

Additional arguments in favor of the modificatio n (|lO| ) 
of the excess free energy (^ will be given in Sect. [V B| . 

Finally, the total free energy 



(11) 



is given by Eqs. (|), {§, §), §), (0), and ( [A211 ). 

Other chemical species can be easily included in our 
analytic model. For example, H2 molecules can be taken 
into account by adding van der Waals a nd HS terms with 
appro priate effective diameters in Eqs. (All), (A16) and 
( [Aig| ) . The generalization of the unbounding term is 



where ric = ne + nu + 7tH2, and wh2 is the average molec- 
ular volume. 



IV. OCCUPATION PROBABILITY FORMALISM 

A. Generalized Saha equation 

The internal free energy (^ can be rewritten in a more 
familiar form. Let us replace the Boltzmann distribution 
Nk, (X Qk, exp{PxK) (which yields the divergent ideal IFF) 
by any real distribution 



where the generalized IFF 



plays role of a normalization constant. Then 



(13) 



(14) 



(15) 



Note that and Z^, in Eqs. (13)-(|l5|) can be multi- 
plied by a common factor. It does not affect the A^^/A^h 
distribution, but should be chosen consistent with the 
ionization equilibrium conditions. 

Although the first sum in Eq. (|lj) is naturally derived 
from the rigorous Eq. (^, it was often omitted in the 
internal free energy ^ or regarded as a part of Fconf [D ■ 
Meanwhile, it has a clear physical meaning: —Nf^lnwK 
is a contribution to the ideal-gas part of the entropy due 
to the correction Wk, to the probability that Kth state is 
occupied. The factors are traditionally called occu- 
pation probabilities, although they do not always have 
direct probability meaning. 

The minimum of the Helmholtz free energy under the 
stoichiometric constraints requires 



OF/dN^ = dF/dNp + dF/dN^. 



(16) 



Fuh = kBTNcinuvii + 71H2WH2) 



(12) 



Separation of i^id from Fconf allows one to rewrite 
Eq. ( p^ ) in the form of the Saha equation: 

n« = npno(ApAo/AH)^w„(gK/2)exp[/3(xK + Adcg)], (17) 

where Adog = is the correction due to the par- 

tial electron degeneracy, and w„ is defined by 

7 rn 1 (^FcQuf dF^Qi^f ^-^conf f-t n\ 

'^^^""'^ = ^ + ^-^- 

This definition is consistent with Eqs. (p^)-(p^, and it 
fixes the above mentioned common factor. These oc- 
cupation probabilities have the same meaning as those 
considered by HM but take into account charged parti- 
cles. Therefore, Eq. (|l|) generalizes Eq. (2.18) of Ref. |] 
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to the case when ionization-recombination processes are 
allowed in the system. 

Equivalently, one may adhere the traditional defini- 
tion, Inwa — —(3dFconi/dNa, and replace tu^, in Eq. (|l^ ) 
by Wk/ {wcWp). However, since proton and electron can- 
not be destroyed by the external fields, we put their "oc- 
cupation probabilities" equal to unity, thus choosing the 
definition ([l8| ) for the atomic occupation factors Wk,. 

An equation equivalent to Eqs. (|l^) and ( p^ ) was de- 
rived by Fontaine et al. (FGVH), who however did 
not introduce the occupation probabilities explicitly. 

If all Wk, were known, then Eq. ( p7| ) would give direct 
solution to the problem. However, since Fconf depends 
on the occupation numbers, Eqs. ( p7| ) and ( p^ ) are to 
be solved together. Nevertheless, this reformulation of 
the problem is useful, because the coupled equations can 
be solved iteratively. First, one chooses an initial value 
of Wn and calculates the particle numbers {Na\ from 
Eq. ([17|) . Then w^, are refined by substituting {N^} into 
Eq. (p^). At low densities, where many excited states are 
populated, this procedure appears to be more efficient 
than alternative schemes |6|, |l3|,^. 

Substituting Eqs. (|ll|) and~(|A2lD into Eq. (|l|) we de- 
compose into five factors. 



(19) 



corresponding to the unbounding of atoms, the Coulomb 
interactions of charged particles, the hard-sphere repul- 
sion, and the corrections due to atom-atom and atom-ion 
attraction. Equation (HG) yields 



„(ub) _ 



exp[-nct;K]- 



(20) 



Expanding Fq and import in powers of particle numbers, 
keeping quadratic terms, and using Eq. (18), we obtain 



Inw^^) 


= -^4^(/3e2)3(ne + np) ( 


1 - 7^7) , 


(21) 


lnzz;(H) 






(22) 


Inw^'") 






(23) 



where a„„' and are the van der Waals constants de- 
fined in Appendix. Our calculations show that the factors 
^y(in) and w^^^ are close to unity, being, thus, unimpor- 
tant. In contrast, w^*^^ is significantly less than unity at 
high densities, even at a relat ively low ionization. 
For the HS repulsion, Eq. ( |A20| ) gives 



(HS) _ (lnz^i°^)(l-W2)-57?^ + 3773 
(1-77)3 

where rj is the filling factor, and 



Inu.?) 



An 



(24) 



(25) 



If 7y <C 1, then w 



(HS) 



W 



(0) 



In practice, the non-linear 



corrections given by Eq. (^4|) may be important. We 
always have it;(^^) < 1, since d^' < d^^,' in our model. 

The formalism can be generalized, e.g., for formation of 
molecules. The dissociation-recombination equilibrium 
is given by 



"H2 = rtH(^H\/2)^Z^2/-^. 



W 7 



(26) 



where Zw2 is the internal molecular partition function, 
generalized through multiplying each /ith term by an oc- 
cupation probability w^^^' (e.g., Eq. (21) of Ref. [^). 
From the equation of chemical equilibrium 



dF/dNm = 2 {dF/dNp + dF/dN^), 



(27) 



we conclude that the molecular occupation probability 
can be defined as 



knTlnw'^^^^ = 2 



fdF 



conf 



dF, 



conf 



OF, 



conf 



(28) 



Keeping the traditional definition of w both for atoms 
and for molecules would result in the same Eq. (pq). 

For the perturbation and HS factors in w^^^ Eqs. 
(p^)-(p5|) remain valid if to use relevant scaling factors; 
the Coulomb factor w^'~"'' should be squared; and the un- 
bounding factor is derived from Eq. 

Inw^^^'"^' = -nH(w^ - wh) - '^H2(wm - WH2) - n-cV^. 

(29) 

As a test example (although marginal to the present dis- 
cussion), we have implemented this approach to molec- 
ular formation, utilizing a simplified treatment of Zi„2 
1^,^. In this case, the unbounding factor ( ^9|) turned 
out to be unimportant, since the degree of dissociation 
is mainly determined by relation between atomic and 
molecular effective HS diameters. 

The situation is very different for the ionization of 
atoms. The repulsion factor ty^^^^ becomes important at 
moderate p and low T, but w^^^^ strongly affects mostly 
excited states at any density (the ground state at high 
density), decreasing rapidly as the conditions for state 
survival become violated. This leads to a physically rea- 
sonable convergence of the IFF at low densities (thus 
radically improving the ionization degree at T > lO"*'^ K) 
and pressure ionization at high densities. 



B. Thermodynamic and optical occupation 
probabilities 

The occupation probability technique described above 
allows one to calculate thermodynamic properties of par- 
tially ionized hydrogen plasma. However, it is still insuf- 
ficient for describing optical properties of the plasma. 

The free energy model presented in Sect. HI and in 



Appendix allows for close configurations of atoms with 
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protons, because the atom-ion repulsion diameters do 
not exceed the quantum-mechanical sizes of atoms. As 
argued in Sect. ||, the close configurations simulate clus- 
ter states (approximately treated as interacting atoms 
and ions in the frames of the chemical picture). This 
approach yields physically plausible EOS. On the other 
hand, specific quantum- mechanical properties (e.g., fre- 
quencies and oscillator strengths of radiative transitions) 
of clusters most likely differ from those of isolated atoms. 
Formation of many different close configurations should 
manifest itself in optics as quasicontinuum. 

Therefore, one should discriminate between the ther- 
modynamic continuum (the states which do not con- 
tribute to the generalized IPF), and the optical con- 
tinuum (the states strongly perturbed by surrounding). 
This dichotomy was first realized by Rogers , who de- 
veloped the concept of the optical and plasma continua 
using the physical picture. A correct account of the qua- 
sicontinuum has been also taken in a recent study of line 
shapes in hydrogen opacities [ p6| . 

The optical continuum can be determined from con- 
sideration of Stark merging of spectral lines of an atom 
affected by plasma microfields This leads to the 

atomic "survival" probabilities w^,, generally different 
from Wk introduced in thermodynamics. Let us call 
i&K the optical occupation probability, to avoid confu- 
sion with w„ . The occupation probabilities based on the 
plasma microiield distribution are in fact the op- 

tical ones. Their implication in thermodynamics leads 
to physically unrealistic EOS (cf. Ref. |1^) due to the 
incorrect treatment of close configurations in the free en- 
ergy. Ionization equilibrium would be equally implausible 
without an ad hoc "pressure ionization" term Indeed, 
substituting the excess free energy (^ into Eq. ([l6| ) we 
would arrive at Eq. (17) with replaced by 



i/;!"^) = exp 



(30) 

Since Vp^ ^ vik, the last (positive) term in Eq. ( |30| ) 
may dominate and yield the "occupation probabilities" 
which grow exponentially with riu- Then one would get 
pressure neutralization instead of pressure ionization at 
high densities |2|] . 

In contrast, Eq. ( [lo| ) leads to the occupation probabil- 
ities, Eq. (p0|), which decrease exponentially with density 
and produce the desired pressure ionization. 

The occupation probabilities given by Eq. (4.71) of 
Ref. |] can be presented (for pure hydrogen plasma) 
in the form ( |30| ) but without the last term. This form 
of Wk cannot be rigorously derived from Eq. (^. Note, 
however, that the leading factor np-Cp^ in the exponent 
occurs due to the Inglis-Teller effect ||2^ which is opti- 
cal but not thermodynamic. Accordingly, these results 
can be used for calculating the fraction of atoms which 
are only slightly perturbed by plasma microfields so that 



they are able to contribute to the atomic opacities. For 
this purpose, we use an expression similar to that in Ref. 
Q . However we take into account that Fpair includes the 
HS term, which implies that the distance between ions 
and atoms cannot be shorter than d^- Thus the HS vol- 
ume should be subtracted from the interaction volume 
V. The latter one has been estimated by several authors 
using different (not always justified) approximations, as 
discussed by HM [||. A reasonable order-of-magnitude 
estimate reads v — |7r(4^„)^. For 7 lowest states this 
estimate is intermediate between more complicated Eqs. 
(4.69) and (4.70) of Ref. g, and for the ground state 
atom it fits the latter exponential with an accuracy of 
9%. Finally, we adopt 



w^ = exp [-f np((4;,)3-d3(r))] 



(31) 



The above considerations emphasize that Wk determine 
solely the optical properties of the plasma and they 
should not be used in the construction of the IPF. 



V. RESULTS AND DISCUSSION 



A. Ionization equilibrium 



The free energy model described in Sect. IH has been 
applied to calculation of the thermodynamic properties 
of plasma using the method of Sect. [V. The ionization 
isotherms are shown in Fig. 1. Light solid curves repre- 
sent the fraction of all H atoms, /h = «h/("-h + "p), 
and dashed lines display the fraction of ground state 
atoms. The results are in general agreement with Ref. 
|l0| but disagree with Ref. If, for example, p = 0.1 g 
cm^'^, Fig. 2 of Ref. |^ shows practically zero ionization 
at T < 10^-^ K, whereas according to Ref. |]l0| there is a 
considerable amount (about 6% by mass) of free protons 
at T = 22 000 K. Our result coincides with the latter one. 
The pressure ionization in our Fig. 1 proceeds smoothly 
at high densities, again in agreement with Ref. ||l0|| , but 
contrary to almost abrupt pressure ionization of Ref. . 

Thick solid and dashed lines in Fig. 1 are obtained 
using Eq. (|3l| ) and show the fraction of those atoms 
whose optical properties are not destroyed by plasma mi- 
crofields, and which therefore should be used in the opac- 
ity calculations. At p ^ 10~^ g cm~^ these curves are in 
good agreement with those in Ref. Ig]. This observation 
suggests a possible explanation to a discrepancy in oc- 
cupation numbers of excited states, recently recognized 
p9[ between OPAL and MDH data: the former ones take 
into account all thermodynamically significant while the 
latter ones only optically identifiable atomic states. 

For comparison, long dashes show the solution of the 
ideal-gas Saha equation including the ground state atoms 
only. These curves reproduce accurately the number 
of ground state atoms at low densities. However, at 
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T > 10"*'^ K the total number density nn can never be de- 
termined in this way, since the excited states become pop- 
ulated and increase the neutral fraction. Above 10^ K, 
the onset of occupation of the excited states produces 
typical "shoulders" on the solid curves. Note that high- 
est populated excited states at low densities are strongly 
affected by microfields and belong to the optical quasi- 
continuum. This explains why the low-density tails of 
the heavy lines lie significantly lower than the light ones. 
This difference is noticeable not only in the pressure ion- 
ization domain p ^ 0.1 g cm"'^ (where the present model 
has a limited applicability), but also at lower densities, if 
the temperature is high enough for population of the ex- 
cited states. On the other hand, when density increases, 
the pressure ionization comes into effect, the neutral frac- 
tion becomes smaller and finally disappears at p > 3 g 
cm^'^. There is a considerable amount of bound species 
at p '--^ 0.1 — 1 g cm~^ (important for thermodynamics) 
which can hardly contribute to the atomic opacities. The 
optical properties of atoms are destroyed at p > 10^^ g 
cm~'^, as is seen from downward bending of the heavy 
lines. 
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FIG. 1. Total neutral fraction (/h) and partial fractions: 
ground state atoms (/hi), optically identifiable bound atoms 
(/opt), and optically identifiable ground state atoms (/opti), 
compared with the ideal gas (Saha) approximation. 



B. Equation of state 

Relative importance of partial contributions to the free 
energy can be estimated by examination of their influence 
on EOS, as illustrated in Fig. 2 for two temperatures. At 
lower temperature, T = 12 600 K, the ionization degree 
is small, and the corrections due to the Coulomb nonide- 
ality of the charged component are practically unimpor- 
tant. The most important corrections are produced by 
the repulsion of atoms (HS) and by the unbounding. The 
corresponding contributions to the pressure are nearly 
equal and become appreciable at p > 0.03 g cm"-^. At 



higher temperature, T = 2 • 10^ K, the ionization is high, 
and the Coulomb nonideality is significant, while the HS 
contribution is nearly negligible. The unbounding, how- 
ever, is important at this temperature also. Note that 
the perturbation (van der Waals) terms are unimportant 
at any T and p shown in Fig. 2. 
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FIG. 2. Partial pressures: ideal gas (id) and configura- 
tional parts due to the Coulomb interactions of the ionized 
fraction (C), unbounding of atoms (ub), strong repulsion at 
short distances (HS), and long-distance attraction (pert). We 
show |Pc| and |Ppcrt| since these parts are negative. 

In Fig. 3 we compare the present EOS (heavy solid 
curves) with those known in the literature. Results of 
FGVH [in are shown by dashes, MDH §] by dot-dashed 
lines, SC Q by light solid lines, and OPAL [|o) by 
dots. The EOS of SC is somewhat softer, while that of 
MDH much stiffer (probably owing to the rapidly increas- 
ing "pressure ionization" term introduced by the latter 
authors). We have terminated our curves at densities 
where the quantum correction in Eq. (||) reached 0.5. At 
higher p the nonideality of the charged component be- 
comes too strong to be treated as perturbation. The 
lowest isotherm corresponds to T = 12 600 K, which is 
below the critical temperature for the plasma phase tran- 
sition reported by SC ||l^ (the density discontinuity seen 
on the corresponding curve). 

Our results have been obtained under the assumption 
that the molecules are completely destroyed. All authors 
agree that it is true at T > 10'* ''' K, however there is 
a great uncertainty concerning the degree of dissocia- 
tion at 4.1 < logioT < 4.5. According to Ref. ||], the 
amount of molecules in this interval is quite insignificant 
at any density, while SC [ p^|JTl] ] and Reinholz et al. p^ ] 
(RRN) found similarly strong but quantitatively different 
recombination at density increasing from 10~^ toward 1 
g cm~^. Quantum molecular dynamics simulations ]3l| ] 
show that there is a significant amount of transient H2- 
like clusters at /cbT = 1 eV and p '--^ 1 g cm~^, although 
their influence on thermodynamic and optical properties 
of plasma is not yet well understood. 
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FIG. 3. Comparison of present EOS with results from the 
literature. 

Figure 4 illustrates the effect of this uncertainty on 
the EOS. The pressure isotherm (T = 15 000 K) cor- 
responding to complete dissociation (heavy solid line) is 
compared with two modified isotherms, obtained through 
replacing in our model the ideal contribution to P by that 
corresponding to the dissociation degree given either by 
SC (dot-dashed line) or by RRN (dashed line). A com- 
parison with the isotherm of SC (drawn by light solid 
line) suggests that the difference between our and SC 
EOS at logj^o ^ ~ 4.5 is mainly caused by the formation 
of dimers, which we neglected. Our testing calculations 
revealed, however, that the dissociation degree in the con- 
sidered p — T domain depends crucially on treatment of 
various molecular excitations as well as on adopted HS 
diameters. Thus the thorny problem of clustering hardly 
can get an unambiguous solution within the chemical pic- 
ture. 

Figure 5 shows the adiabatic temperature gradient 
Vad = {d\nT/d\nP)s, a significant quantity sensitive to 
details of the free energy model. Thermal radiation starts 
to dominate in thermodynamics for the hottest isotherm 
at low density, causing the decrease of Vad. Other de- 
pressions of Vad are explained by the increase of the spe- 
cific heat in the regions of partial ionization, where the 
internal energy is affected by the strongly T-dependent 
ionization degree. Our data (heavy curves) are com- 
pared with the tables of FGVH |l| (dashes), SC Q 
(dot-dashed lines, "table") and OPAL g] (dots). Light 



solid lines ("formula") are obtained by substituting the 
SC tabulated quantities P, S, T, a = (9 log p/9 log T)p 
and b — {dlogS/d\ogT)p into the thermodynamic iden- 
tity Vad = -aPV/{bST). We have plotted only the 
latter lines for the two highest isotherms, for which the 
"table" and "formula" values are in a good agreement. 
The discrepancy between these values in the low-T high- 
p domain indicates the lack of thermodynamic consis- 
tency |ll] caused by a high sensitivity of the complicated 
numerical approach to accidental small errors in mini- 
mization and differentiation procedures. Although our 
simplified analytic model is less scrupulous in details, its 
advantage is that it is free of such inconsistencies. 
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FIG. 4. Influence of the uncertainty in the dissociation 
degree at high density on the EOS. 

The most crucial test for the validity of our model may 
be provided by a comparison with the advanced results 
based on the physical picture and employed in the OPAL 
opacity library ||30| . They do not cover the most interest- 
ing region where other models reveal major discrepancies, 
but the available data are in satisfactory agreement with 
our model. 

Sequences of shallow depressions in the isotherms T > 
10^ K at p > 10~^ g cm~^ indicate successive pressure 
destruction of excited atomic states. Their physical real- 
ity remains an open question. Note that the SC results 
reveal analogous oscillatory behavior (which is probably 
inherent to the models of such a type), which however is 
not observed in the OPAL data. 
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introduce the optical occupation probabilities which de- 
termine the neutral fraction visible in atomic lines. The 
non-perturbed atomic fraction given by these probabili- 
ties agrees with that of Ref. |^ at the densities available 
in laboratory. 

The equation of state obtained from our model is com- 
pared with the results of other authors. The best agree- 
ment is achieved with the OPAL data (in the p — T 
region where they are available). Since the equation of 
state employed in OPAL is based on the physical picture 
of plasma, completely different from our model, we regard 
this agreement as an indirect confirmation of the validity 
of our approach. Its generalization to higher densities 
and higher atomic numbers is being performed. 
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FIG. 5. Adiabatic temperature gradient compared with 
the data of FGVH, SC, and OPAL. The curves are labeled by 

VI. SUMMARY 

We have developed an analytic free energy model for 
partially ionized hydrogen plasma in the framework of 
the chemical picture. The model describes thermody- 
namic properties of the plasma at T > 10^ K and 
p ^ 0.1 — 1 g cm^'^. In particular, it can be used in the 
studies of DA white dwarf and neutron star atmospheres. 

The occupation probability formalism, first introduced 
by Fermi Ig^l and further developed by Hummer and Mi- 
halas is generalized to take proper account of the 
effects of partial ionization. Free energy minimization is 
obtained by a generalized Saha equation which is solved 
by an iterative algorithm involving the modified occu- 
pation probabilities. Calculated ionization degree differs 
from that obtained previously in Ref. but qualita- 
tively agrees with the results of Ref. jl^ . We argue that 
the relatively high neutral fraction obtained in ReL pc| ] 
and in our present work at p '--^ 0.1 — 1 g cm~^ cannot be 
pronounced in atomic opacities, since it takes cumulative 
account of atomic and cluster states, the latter ones con- 
tributing to the optical quasicontinuum (Sect. [TV B| ). We 
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APPENDIX: FREE ENERGY MODEL BASED ON 
INTERPARTICLE PAIR POTENTIALS 

Our model correction to the excess free energy Fpair 
due to binary interactions among atoms and protons is 
mainly based on Refs. 1^,^. Some modifications, how- 
ever, are introduced to improve the physical consistency 
of the model, and analytic fits are presented. 

The total excess free energy consists of the reference 
part FhS: treated in the hard-sphere approximation, and 
perturbation parts Fpcrt required to include the attrac- 
tive (van der Waals) interactions. 

1. Atom-atom interactions 

We treat interactions among neutral species using the 
fluid perturbation theory of Weeks, Chandler, and An- 
dersen (WCA). An effective potential ^(r) is sepa- 
rated into the reference and perturbation parts, 

0(r) = (?!)rcf (r-) + '?!)pcrt(?-), (Al) 

where 0rcf is a 'purely repulsive finite-range potential 
which acts at distances r < r*, where is the mini- 
mum point of (/)(r). Accordingly, the free energy splits 
into the reference and perturbation parts. First we con- 
sider the ground state atoms and adopt the interatomic 
potential 
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(A2) 



where 0g^^ (r) is given by Eq. (2) of Rcf . @, and ^iorU^) 
is a correction which acts at short distances r < 1.5 A. 
The correction is required to make (/)(r) go to infinity at 
r ^ 0. Unhke SC ||l^ who introduced HS cores with 
fixed ad hoc diameters, we use the scaled mean electric 
potential of the ground state H atom: 



.H,M4(l + I)e-/-. 



(A3) 



This choice seems reasonable since it yields the Coulomb 
repulsion at r — > 0. We use the effective He screening 



length 



Ob, where ob is the Bohr radius, for 



the ground state atom. 

The perturbation part of the free energy is 



2 



0pcrt('')5rcf(»')d^f, 



(A4) 



where ^rcf is the pair correlation function of the refer- 
ence system. For the HS reference system at low densi- 
ties, grefii") can be replaced by at r < d, and 1 at r > d, 
where d is the HS diameter. Then 



(A5) 



where co/Zcb = 20.2 K and i?o — 8.6aB for the ground 
state atom. 

The diameter d can be determined from the WCA self- 
consistency condition, which involves radial distribution 
function. However, according to Ref. ||34[|, the WCA 
value of d at low p is close to the Barker's |35[ value 



[l-exp(-/50,rf(r))]dr. 



(A6) 



We adopt d = d-Q aX p <^ 1 g cm^'^, and propose the 
expression d ^ d^ cxp{—p/2 g cm^"^) at p ^ 1 g cm^'^. 
The latter expression fits exact numerical results Q with 
an error of about 5%. 



We have used the reference part of the potential ( [A2| ) 
in Eq. (|A6|) and fitted the result by the formula 



do 



1 + ln(l + ci\/i) + 



t 



1 + ln(l + C2<) 



(A7) 



where do = Gqb, ci 4, = 0.5, and i = r/(3. 25-10'' K) 
for the ground state atom. 

We use simple scaling of the potentials for excited 
states. For two atoms with principal quantum numbers 
n and n', we write 



'V«'/'H(?'/Sn'n), 



(A8) 



where s„'„ ~ (/„ + ln')/i2h) is the scali ng f actor, and Z„ 
is an average atomic size. The scaling ( [A8| ) ensures the 
correct Coulomb repulsion at short distances. We esti- 
mate l'^ as the quantum-mechanical expectation value of 



yields 



averaged over the quantum numbers (I, m), which 



= Ob n V(7n2-|-5)/4 



(A9) 



The scaling does not reduce to just multiplying djj by s, 
but implies simultaneous scaling of the temperature: 

(T) = T). (AlO) 

The perturbation terms for different states are additive: 

i^(H) _ _ AT AT . I\f (All) 



where the van der Wa als constants a^^.' are determined 
by the scaled Eq. (A5), 



= -^L' £0 Y {Rl - dU^^.'T)) . (A12) 



2. Ion-atom interactions 

Following SC we describe the polarization inter- 
action outside the core by the screened dipolar potential: 



II + r2 



(A13) 



where a„ is the average polarizability of an atom with 
the principal quantum number n, and td is the scre ening 
length. The rms size In in the denominator of Eq. ( A13| ) 
is intermediate between two different values of the polar- 
ization radius used in Ref. ||l0|| . Furthermore, the polar- 
izability of H atom in the state (nlm) averaged over (Im) 
numbers [|3q| can be fitted, with an error of 2%, by 



0.85^; 



(AM) 



At short separations, the Coulomb re puls ion should pre- 
vail. Therefore, in analogy with Eqs. ( A2), we adopt the 
interaction potential 



(A15) 



where 0corc is given by Eq. (A3) with a = ^SnOB and the 
scaling factor s„ = In/ h- 

The potential (A15) is separate d th en into the refer- 
ence and perturbation parts, Eq. ( |Al[) , according to the 
WCA prescription. The pe rtur bation free energy is cal- 
culated analogously to Eq. (A4) (with rip instead of 
on the right-hand side). This gives 



where 



2/2 



(A16) 



(A17) 
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Here e^Z^w* is the integral of 0p"j?j(r) over space, 
{—e'^/li^)4>!t is the minimum of (/lin, and d^iT) = 
lKd^{sK,T) is the HS diameter determined by Eq. ( [A^ ) 
(nearly independent of td). We have obtained the fits 



4- 1.7.T 
1 + 1.7a;' 



4>*{x) 



0.155- 0.0212j;^ 
1 + 0.34x2 



(A18) 



and d*(r) is given by Eq. (A?) with do = 0.615, ci — 
0.71, C2 = 0.75, and t = r/(2.15 • 10^ K). 



3. Hard sphere contribution 

The HS diameters depend on atomic states and dif- 
fer for interactions with atoms and ions. Thus we have 
a non-additive HS mixture. Such mixtures can be de- 
scribed by the van der Waals one-fluid model, which is 
reasonably accurate for effective filling factors 77 < 0.3 
|p7| . In the spirit of this model, we define 



and use the Carnahan-Starling [Q formula 



(A19) 



(A20) 



Finally, the excess free ene rgy a ss ociat ed wit h th e pair 
potentials is given by Eqs. ( [A20| ), (^nj) and ( |aT^ ): 



Fm + F^ 



(H) 



pert 



F, 



(in) 



pert • 



(A21) 
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